arXiv: 1506.05121 v2 [astro-ph.GA] 24Jul2015 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 27 July 2015 (MN I4 TbX style file v2.2) 


Radiation pressure driving of a dusty atmosphere 

Benny T.-H. Tsang and Milos Milosavljevic 

Department of Astronomy, University of Texas at Austin, Austin, TX 78712, USA 


27 July 2015 


ABSTRACT 

Radiation pressure can be dynamically important in star-forming environments such as ultra- 
luminous infrared and submillimeter galaxies. Whether and how radiation drives turbulence 
and bulk outflows in star formation sites is still unclear. The uncertainty in part reflects the 
limitations of direct numerical schemes that are currently used to simulate radiation trans¬ 
fer and radiation-gas coupling. An idealized setup in which radiation is introduced at the 
base of a dusty atmosphere in a gravitational field has recently become the standard test for 
radiation-hydrodynamics methods in the context of star formation. To a series of treatments 
featuring the flux-limited-diffusion approximation as well as a short-characteristics tracing 
and Ml closure for the variable Eddington tensor approximation, we here add another treat¬ 
ment that is based on the Implicit Monte Carlo radiation transfer scheme. Consistent with 
all previous treatments, the atmosphere undergoes Rayleigh-Taylor instability and readjusts 
to a near-Eddington-limited state. We detect late-time net acceleration in which the turbulent 
velocity dispersion matches that reported previously with the short-characteristics-based ra¬ 
diation transport closure, the most accurate of the three preceding treatments. Our technical 
result demonstrates the importance of accurate radiation transfer in simulations of radiative 
feedback. 

Key words: star: formation - ISM: kinematics and dynamics - galaxies: star formation - 
radiative transfer - hydrodynamics - methods: numerical 


1 INTRODUCTION 


The forcing of gas by stellar and dust-reprocessed radiation has 
been suggested to reduce star formation efficiency and drive super¬ 


sonic turbulence and large-scale outflows in galaxies (e.g., Thomp 


et al.|2005 2015||Murray et al.|2010[|2011||Faucher-Giguere 

et al.|2013[|Kuiper et al.|2015] . Generally the net effect of radiation 
pressure is to counter the gravitational force and modulate the rate 
of infall and accretion onto star-forming sites. In its most extreme 
presentation, radiation pressure accelerates gas against gravity so 
intensely that the gas becomes unbound. For example, |Geach et al.| 
( [2014] ) have recently suggested that stellar radiation pressure drives 
the high-velocity, extended molecular outflow seen in a starburst 
galaxy at z = 0.7. Theoretical and observational evidence thus 
suggests that radiation may profoundly influence the formation and 
evolution of star clusters and galaxies. While direct radiation pres¬ 
sure from massive young stars may itself be important in some, 
especially dust-poor environments (Wise et al.|2012| l, the trapping 
of the stellar radiation that has been reprocessed by dust grains into 
the infrared (IR) should be the salient process enabling radiation 
pressure feedback in systems with the highest star formation rate 
densities. 


Usually, the amplitude of radiative driving of the interstellar 
medium (ISM) in star-forming galaxies is quantified with the av¬ 
erage Eddington ratio defined as the stellar UV (or, alternatively, 
emerging IR) luminosity divided by the Eddington-limited lumi¬ 


nosity computed with respect to the dust opacity. However, in re¬ 
ality, the ISM is turbulent and dust column densities vary widely 
between different directions in which radiation can escape. The lo¬ 
cal Eddington ratio along a particular, low-column-density direc¬ 
tion can exceed unity even when the average ratio is below unity. 
[Thompson & Krumholz|(2014] argue that this is sufficient for radi¬ 
ation pressure to accelerate gas to galactic escape velocities. |An-| 
|drews & Thompson| (201 1| > surveyed star forming systems on a 
large range of luminosity scales, from star clusters to starbursts, 
and found that their dust Eddington ratios are consistent with the 
assumption that radiation pressure regulates star formation. 

|Murray et al.| (2005] highlighted the importance of radiation 
momentum deposition on dust grains in starburst galaxies. They 
showed that the Faber-Jackson relation L oc a 4 and the black hole 
mass-stellar velocity dispersion relation Mbh oc a 4 could both be 
manifestations of self-regulation by radiation pressure. Thompson 
|et al.| ( [20Q5] l argued that radiation pressure on dust grains can pro¬ 
vide vertical support against gravity in disks of starburst galaxies 
if the disks are optically thick to the reprocessed IR radiation. In 
a study of giant molecular cloud (GMCs) disruption, |Murray et al.| 
(2010] found that radiation pressure actually dominates in rapidly 
star-forming galaxies such as ULIRGs and submillimeter galaxies. 
|Hopkins et al.|(2010] identified a common maximum stellar mass 
surface density £ max ~ 10 11 Mq kpc~ 2 in a variety of stellar sys¬ 
tems ranging from globular clusters to massive star clusters in star- 
burst galaxies and further to dwarf and giant ellipticals. These sys- 
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terns spanned ~ 7 orders of magnitudes in stellar mass and ~ 5 
orders of magnitude in effective radius. The universality of maxi¬ 
mum stellar mass surface density can be interpreted as circumstan¬ 
tial evidence for the inhibition of gaseous gravitational collapse by 
radiation pressure. 

The preceding studies were based on one-dimensional or oth¬ 
erwise idealized models. To understand the dynamical effects of 
radiation pressure in a dusty ISM. however, multi-dimensional ra¬ 
diation hydrodynamics (RHD) simulations are required. One spe¬ 
cific setup has emerged as the testbed for radiation hydrodynam¬ 
ics numerical methods used in simulating the dusty ISM, specifi¬ 
cally in the regime in which the gas (assumed to be thermally cou¬ 
pled to dust) is approximately isothermal and susceptible to com¬ 
pressive, high-mach-number perturbations. |Krumholz & Thomp-| 
|son| (2012| hereafter KT12) and |Krumholz & Thompson (20 1 3[ 
hereafter KT13) designed a two-dimensional model setup to inves¬ 
tigate the efficiency of momentum transfer from trapped IR radi¬ 
ation to a dusty atmosphere in a vertical gravitational field. Using 
the flux-limited diffusion (FLD) approximation in the ORION code 
(Krumholz et al.|2007| l, they found that the optically thick gas layer 
quickly developed thin filaments via the radiative Rayleigh-Taylor 
instability (RTI). The instability produced clumping that allowed 
radiation to escape through low-density channels. This significantly 
reduced net momentum transfer from the escaping radiation to the 
gas, and the gas collapsed under gravity at the base of the compu¬ 
tational box where radiation was being injected. 

|Davis et al.| < |20l4] hereafter D14) then followed up by sim¬ 
ulating the same setup with the ATHENA code (Davis et al.|20 1 2\ 
using the more accurate variable Eddington tensor (VET) approx¬ 
imation. They constructed the local Eddington tensor by solving 
the time-independent radiative transfer equation on a discrete set 
of short characteristics (Davis et al.|2012} . Similar to the simula¬ 
tions of KT12, those of D14 developed filamentary structures that 
reduced radiation-gas momentum coupling. However, in the long¬ 
term evolution of the radiation-pressure-forced atmosphere, D14 
detected significant differences, namely, the gas continued to accel¬ 
erate upward, whereas in KT12, it had settled in a turbulent steady 
state confined near the base of the box. D14 interpreted this dif¬ 
ference of outcome by referring to an inaccurate modeling of the 
radiation flux in the optically thick-to-thin transition in FLD. 

|Rosdahl & Teyssierj (2015| hereafter RTI5), simulated the 
same setup with the new RAMSES-RT RHD code using the com¬ 
putationally efficient Ml closure for the Eddington tensor. This 
method separately transports radiation energy density and flux and 
assumes that the angular distribution of the radiation intensity is a 
Lorentz-boosted Planck specific intensity. One expects this to pro¬ 
vide a significant improvement of accuracy over FLD, however still 
not approaching the superior accuracy of the short characteristics 
closure. The Ml results are qualitatively closer to those obtained 
with the FLD than with the short-characteristics VET. RTI5 ar¬ 
gue that the differences between FLD and Ml on one hand and the 
short characteristics VET on the other may be more subtle than sim¬ 
ply arising from incorrectly approximating the flux at the optically 
thick-to-thin transition. 

In this paper, we revisit the problem of radiative forcing of 
a dusty atmosphere and attempt to reproduce the simulations of 
KT12, D14, and RT15, but now with an entirely different numerical 
scheme, the implicit Monte Carlo (IMC) method of Abdikamalov| 
|et al.| f2012) originally introduced by Fleck & Cummings t'1971| >. 
The paper is organized as follows. In Section[2]we review the equa¬ 
tions of radiation hydrodynamics and the IMC method. In Section[3] 
we then assess the reliability of our approach in a suite of standard 


radiation hydrodynamics test problems. In Section [4] we describe 
the simulation setup and details of numerical implementation. We 
present our results in Section[5]and provide concluding reflections 
in Section[6] 


2 CONSERVATION LAWS AND NUMERICAL SCHEME 


We start from the equations of non-relativistic radiation hydrody¬ 
namics. The hydrodynamic conservation laws are 


! + V.( P v),0, 


0 ) 


+ V ■ (pvv) + VP = pg + S, 


<9pv 
dt 

^+V-[(pP + P)v] = pvg + cS 0 , 

where p, v, and P are the gas density, velocity, and pressure of the 
gas and 


( 2 ) 


( 3 ) 


E = e + 


( 4 ) 


is the specific total gas energy defined as the sum of the specific 
internal and kinetic energies of the gas. On the right hand side, g is 
the gravitational acceleration and S and cSa are the gas momentum 
and energy source terms arising from the coupling with radiation. 

The source terms, written here in the lab frame, generally de¬ 
pend on the gas density, thermodynamic state, and velocity. We 
here investigate a system that steers clear of the dynamic diffusion 
regime, namely, in our simulations tv/c 1 is satisfied at all 
times. Here, r < 10 2 is the maximum optical depth across the box 
and the velocity is non-relativistic v/c < 10 -4 . Therefore, we can 
safely drop all 0(v/c ) terms contributing to the momentum source 
density and can approximate the lab-frame momentum source den¬ 
sity with a velocity-independent gas-frame expression 


S = 


c/ dt / [ fc ( £ ) J ( e ’ n ) - i( £ > n )] 1 


( 5 ) 


Here, J(e, n) is the specific radiation intensity, e is the photon en¬ 
ergy, n is the radiation propagation direction, c is the speed of 
light, dfl is the differential solid angle in direction n, and k(e) and 
j (e, n) are the total radiation absorption and emission coefficients. 
The coefficients are also functions of the gas density p and temper¬ 
ature T but we omit these parameters for compactness of notation. 

Our energy source term includes the mechanical work per unit 
volume and time v ■ S that is performed by radiation on gas. Since 
the gas-frame radiation force density is used to approximate the 
lab-frame value, the source term is correct only to 0(v/c) 


cS 0 = 


I?<! 


d,VL [k(e)I(e, n) - j(e, n)] + v ■ S. 


(6) 


Because in this scheme radiation exerts force on gas but gas does 
not on radiation, the scheme does not conserve energy and momen¬ 
tum exactly. However, it should be accurate in the non-relativistic, 
static-diffusion limit; we test this accuracy in Section [3] We split 
the absorption and emission coefficients by the nature of radiative 
process 


fc(e) = fca(e) + fcs(e), 


( 7 ) 


j(e,n) = j a (e) + j s (e,n). (8) 

The subscript ‘a’ refers to thermal absorption and emission, and ‘s’ 
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refers to physical scattering (to be distinguished from the effective 
scattering that will be introduced in the implicit scheme). 

Equations CHD couple to the radiation subsystem via the radi¬ 
ation source terms defined in Equation and ID Assuming local 
thermodynamic equilibrium (LTE), the radiation transfer equation 
can be written as 

1 d/(e n) + n V/(e, n) = fc a (e)S(e) - k(e)I(e, n) 
c at 

+ js(e, n) + j ex t(e,n) (9) 

where 13(e) is the Planck function at temperature T and j ex t(e, n) 
is the emissivity of external radiation sources. Note that the j s term 
depends implicitly on the specific intensity I which makes Equa¬ 
tion ID an integro-differential equation. 

Since the absorption and emission coefficients depend on the 
gas temperature, and the temperature in turn evolves with the ab¬ 
sorption and emission of radiation, the system is nonlinear. We 
solve the system by operator-splitting (Section ED’ by replacing 
a portion of absorption and emission with effective scattering (thus 
making the solution implicit; Section [2~2} , and by discretizing the 
radiation field with a Monte-Carlo (MC) scheme (Section[2.3|l. 


2.1 Operator-splitting scheme 

Our numerical method is based on the adaptive-mesh refinement 
(AMR) code FLASH ( |Fryxell et al.|200()||Dubey et al.|2(X)8V ver¬ 
sion 4.2.2. We use operator-splitting to solve Equations HD and 
ID in two steps: 

(i) Hydrodynamic update: Equations HD without the radiation 


source terms 



(10) 

ppp + V ■ (pvv) + VP = pg, 

(ID 


dpE 

dt 


+ V ■ \{pE + P) v] = pv ■ g 


( 12 ) 


are solved using the HYDRO module in FLASH. 

(ii) Radiative transport and source deposition update: Equation 
0 coupled to the radiative momentum 



(13) 


and energy 



(14) 


deposition equations is solved with the implicit method that we pro¬ 
ceed to discuss. 


2.2 Implicit radiation transport 

Under LTE conditions, the tight coupling between radiation and gas 
is stiff and prone to numerical instability. This limits the applicabil¬ 
ity of traditional, explicit methods unless very small time steps are 
adopted. The method of Fleck & Cummings (1971) for nonlinear 
radiation transport relaxes the limitation on the time step by treating 
the radiation-gas coupling semi-implicitly. Effectively, this method 
replaces a portion of absorption and immediate re-emission with 
elastic scattering, thus reducing the amount of zero-sum (quasi¬ 
equilibrium) energy exchange between gas and radiation. Numer¬ 
ous works have been devoted to investigating the semi-implicit 


scheme’s numerical properties. |Wollaber| ( |2008| provides a detailed 
description of the approximations made and presents a blueprint for 
implementation and stability analysis. Cheatham 2010 1 provides 
an analysis of the truncation error. Recently, Roth & Kasen|f20T5) > 
developed variance-reduction estimators for the radiation source 
terms in IMC simulations. 

In this section, we describe a choice of formalism for solving 
the coupled radiative transport and source term deposition equa¬ 
tions. The radiation transport equation is revised to reduce thermal 
coupling between radiation and gas by replacing it with a pseudo¬ 
scattering process. The detailed derivation of the scheme can be 
found in |Fleck& Cummings| l |197l) and |Abdikamalovet al.| ( |2012^ . 
Here we reproduce only the main steps. Our presentation follows 
closely|Abdikamalov et al.|((2012]) and the approximations are as in 
|Wollaber|d2008| >. 


Given an initial specific intensity J(e, n, t n ) and gas specific 
internal energy e(t n ) at the beginning of a hydrodynamic time step 
t n , our goal is to solve Equations ID and ( |14| ) to compute the time- 
advanced values J(e, n, t n+1 ) and e(t n+i ) at the end of the time 
step t n+1 = t n + At, where At is the hydrodynamic time step. 
During this partial update, we assume that the gas density p and 
velocity v remain constant. For mathematical convenience we in¬ 
troduce auxiliary parametrizations of the thermodynamic variables: 
the gas internal energy density 


Ug = pe, 


(15) 


the energy density that radiation would have if it were in thermo¬ 
dynamic equilibrium with gas 


47t f 00 

Mr = — / B(e)dt (16) 

c Jo 

(for compactness of notation and at no risk of confusion, we do not 
explicitly carry dependence on the gas temperature T), the normal¬ 
ized Planck function 


Me) = B < c > _ 

U Trr f 0 °°B(e)de’ 

the Planck mean absorption coefficient 


(17) 


, _ f 0 °° k a {e)B(e)de 

P fcTBf)de 


(18) 


and a dimensionless factor, /3, quantifying the nonlinearity of the 
gas-radiation coupling 



(19) 


At a risk of repetition, we emphasize that the M r defined in Equation 
( | 1 6| > is not the energy density of the radiation field; it is simply as 
an alternate parametrization of the gas internal energy density. The 
absorption coefficient and the gas-temperature Planck function, in 
particular, can now be treated as functions of u r . 

Taking the physical scattering to be elastic, Equations (|9]l and 
( | 1 3| ), and |l4j can be rewritten as 

1 dl(e n) + R VI(e, n) = fe a (e)6(e)cw r - fc a (e)/(e, n) 

c at 

- k s (e)I(e, n) + j s (e, n) 

+ jext(e, n), (20) 

fldt + ckpUl = J de j d ^Me)/(e, n ). (21) 
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We linearize the equations in u r and 7(e, n) and denote the corre¬ 
sponding constant coefficients with tildes, 


1 dl(t, n) 


= -n ■ V/(e, n) + fe a (e)&(e)cu r - fc a (e)/(e, n) 
- k B (e)I(e,n) + j s (e, n) + j ext (e. n), (22) 


1 du T 
J dt 

The scattering emission coefficient is 


d0fc a (e)7(e, n). 


-f*/ 

fission coefficii 
j s {e, n) = J dfi' fc s (e)S(e,n,n')/(e, n'), 


(23) 


(24) 


where S(e, n, n') is the elastic scattering kernel. We evaluate the 
constant coefficients explicitly at t n , the beginning of the time step. 
Next, for t n ^ t ^ t n+1 , we expand u r to the first order in 

time 


u T (t) ~ u" + (f — 


(25) 


where u" = u T (t n ) and u' T n = du T /dt (t n ). It is worth noting 
that the implicitness in TMC’ refers to that introduced in Equation 
l |25| . Substituting u T (t) from Equation {25} into {23}, solving for 
w'", and substituting the result back into Equation ({25}, we obtain 


Mr = fu? + 


1 -f 

cfc D 


f*/ 


dfl fc a (e)/(e, n) 


where / is a time-dependent factor 

1 


/(*) = 


1 + (t — f n )/3cfc p 


(26) 


(27) 


Since it is desirable to work with time-independent coefficients, 
we approximate /(f) with the so-called Fleck factor that remains 
constant during the time step 

1 


/- 


1 + aAtfickp 


(28) 


where 0 ^ a ^ 1 is a coefficient that interpolates between the 
fully-explicit (a = 0) and fully-implicit (a = 1) scheme for updat¬ 
ing u r . For intermediate values of a, the scheme is semi-implicit. 
The scheme is stable when 0.5 ^ a ^ 1 (Wollaber|2008) . 

We substitute u r from Equation {26} into Equation {20} to ob¬ 
tain an equation for I (e, n) in the form known as the implicit radi¬ 
ation transport equation 


1 dl(e, n) 
c dt 


+ n ■ V7(e, n) = 


fc ea (e)6(e)cM™ - fe ea (e)7(e, n) - k eB (e)I(e, n) 

+ Me)Ke) r d , r d ^,~ k ') 

k P Jo J 

- fe s (e)7(e,n) +/ s (e,n) +/ext(e,n), 


(29) 


where the effective absorption and scattering coefficients are 
fcea(e) =/fea(e), (30) 

fces(e) = (1 - /) fca(e). (31) 

Equation {29} admits an instructive physical interpretation. 
The first two terms on the right-hand side represent the emission 
and absorption of thermal radiation. Direct comparison with Equa¬ 
tion {20} shows that both terms are now a factor of / smaller. 
The following two terms containing k eB are new; their functional 
form mimics the absorption and immediate re-emission describing 


a scattering process. Meanwhile, the physical scattering and exter¬ 
nal source terms have remained unmodified. 

Since the effective absorption fc ea and scattering k eB coeffi¬ 
cients sum to the actual total absorption coefficient fc a , we can in¬ 
terpret Equation {29} as replacing a fraction (1 — /) of absorp¬ 
tion and the corresponding, energy-conserving fraction of emission 
by an elastic psedo-scattering process. The mathematical form of 
the Fleck factor can be rearranged to make this physical interpreta¬ 
tion manifest. Assuming an ideal gas equation of state, the radiative 
cooling time is 


tcool - 


4 

cftkp 


(32) 


and the Fleck factor is 


1 + 4aAf/f coo i' 


(33) 


When Af/f coo i S> 1 so that / <C 1, the absorbed radiation is re¬ 
radiated within the same time step at zero net change in the gas 
energy density; the only change is randomization of the radiation 
propagation direction. The stability of the scheme rests precisely on 
this reduction of the stiff thermal coupling. However, excessively 
large time steps can still produce unphysical solutions jWollaber] 

|2008l >. 

After the radiation transport equation has been solved using 
the IMC method (see Section [23}, the net momentum and energy 
exchange collected during the radiative transport solve, which read 


S = 


J dek(e) J dQ J dtl(e, n) i 


and 


cS 0 

= — 4:TTCU" 


1 / 
At J 0 


+v ■ S, 


+ ~Kt f 0 d£ ke3,{e) J dQ. dtl(e.n) 


(34) 


(35) 


are deposited in the hydrodynamic variables. Therefore step (ii) in 
our operator splitting scheme has now been further split into two 
sub-steps: 

(ii / ) Radiative transport and hydrodynamical source term col¬ 
lection '.: Solve Equation {29} with the IMC method while accumu¬ 
lating the contribution of radiative processes to gas source terms as 
in Equations {34} and {35}. 

(ii”) Hydrodynamical source term deposition: Update gas mo¬ 
mentum and energy density using Equations {13} and {14}. 


2.3 Monte Carlo solution 

The transition layer between the optically thick and thin regimes 
strains the adequacy of numerical radiation transfer methods based 
on low-order closures. In this transition layer, the MC radia¬ 
tive transfer method should perform better than computationally- 
efficient schemes that discretize low-order angular moments of 
Equation {9}. In the MC approach, one obtains solutions of the 
radiation transport equation by representing the radiation field 
with photon packets and modeling absorption and emission with 
stochastic events localized in space and/or time. This permits accu¬ 
rate and straightforward handling of complicated geometries and, 
in greater generality than we need here, angle-dependent physi¬ 
cal processes such as anisotropic scattering. The specific intensity 
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7(e, n) is represented with an ensemble of a sufficiently large num¬ 
ber of MCPs0 

In the radiation transfer update, starting with the radiation field 
at an initial time t n , we wish to compute the coupled radiation-gas 
system at the advanced time t n+1 = t n + At. Our MC scheme 
follows closely that of Abdika malov et aT] (2012| > and Wollaber 
(20081. The radiation field is discretized using a large number of 
Monte Carlo particles (MCPs), each representing a collection of 
photons. We adopt the grey approximation in which we track only 
the position and the collective momentum of the photons in each 
MCP. MCPs are created, destroyed, or their properties are modified 
as needed to model emission, absorption, scattering, and propaga¬ 
tion of radiation. 

If a finite-volume method is used to solve the gas conserva¬ 
tion laws, the physical system is spatially decomposed into a finite 
number of cells. For the purpose of radiative transport, gas prop¬ 
erties are assumed to be constant within each cell. Particles are 
created using cell-specific emissivities. Each MCP is propagated 
along a piecewise linear trajectory on which the gas properties (ab¬ 
sorption and scattering coefficients) are evaluated locally. Our MC 
scheme computes an approximation to the solution of Equation j29j 
in two steps: by first creating MCPs based on the emissivities and 
boundary conditions, and then transporting MCPs through space 
and time. 


2.3.1 Thermal emission 

In Equation \29\ , the term k ea bcu " on the right-hand side is the 
frequency-dependent thermal emissivity. Assuming that thermal 
emission is isotropic (fc ea is angle-independent), the total thermal 
radiation energy emitted by a single cell of gas AS can be calcu¬ 
lated as 

POO 

AS = AtvAtAV / k ea (e)B(e)de, (36) 

Jo 

where At is the time step size, AV is the cell volume, and B{f) is 
the Planck function at the gas temperature T(t n ). Since we further 
assume that the opacity is grey (independent of e) then 

AS = cAtAV~k es .u?. (37) 

The net momentum exchange due to thermal emission is zero be¬ 
cause the thermal radiation source is isotropic. 

We specify that in thermal emission, N new MCPs are cre¬ 
ated in each cell in each time step. The energy carried by each new 
MCP is then AS /N. The emission time of each such MCP is sam¬ 
pled uniformly within the interval [t n ,t n+1 ]. The spatial position 
of the MCP is sampled uniformly within the cell volume and the 
propagation direction is sampled uniformly on a unit sphere. Ev¬ 
ery MCP keeps track of its time ti, current position r,., momentum 
Pi, and fraction of the energy remaining since initial emission q, 


1 One disadvantage of the MC scheme is low computational efficiency in 
the optically thick regime where the photon mean free path is short. Effi¬ 
ciency in such regions can be improved by applying the diffusion approxi¬ 


mation (Fleck & Canfield 1984||Gentile|200Tj |~Densmore et a l.|2007| l. Re- 


cently, [Abdikamalov et al. (2012) interfaced the IMC scheme at low optical 
depths with the Discrete Diffusion Monte Carlo (DDMC) method of |Dens- | 
more et ak| (2007) at high optical depths. This hybrid algorithm has been 
extended to Lagrangian meshes (Wollaeger et al. 2013 1 . In the present ap¬ 
plication, the optical depths are relatively low and shortness of the mean 
free path is not a limitation. 


where the index i ranges over all the MCPs active in a given hydro- 
dynamic time step. Newly created MCPs are added to the pool of 
MCPs carried over from previous hydrodynamic time steps. 


2.3.2 Absorption 

To minimize noise, we treat absorption deterministically. This ‘con¬ 
tinuous absorption’ method is a variance reduction technique com¬ 
mon in practical implementations of IMC (Abdikamalov et ak| 
|2012||Hykes & Densmore|2009) . Specifically, when an MCP trav¬ 
els a distance cSti inside a cell with absorption coefficient fc ea , its 
momentum is attenuated according to 

Pi (ft + 5ti) = p i(ti)e -SeaC,5ti , (38) 

where we denote an arbitrary time interval with 5ti to distinguish 
it from the hydrodynamic time step At. 


2.3.3 Transport 

In a single hydrodynamic time step, the simulation transports the 
MCPs through multiple cells. The MCP-specific time remaining 
until the end of the hydrodynamic time step is t n+1 — ti. For each 
MCP, we calculate or sample the following four distances: 


(i) The free streaming distance to the end of the hydrodynamic 
time step dt = c(t n+1 — ti). 

(ii) Distance to the next scattering event assuming cell-local 
scattering coefficients 


d s = — - 


In? 


k s + ke. 


(39) 


where £ is a random deviate uniformly distributed on the interval 

( 0 , 1 ]- 

(iii) Distance to near-complete absorption d a defined as the dis¬ 
tance over which only a small fraction ? m in = 1CP 5 of the initial 
energy remains. 

(iv) Distance to the current host cell boundary db. 


We repeatedly update the four distances, select the shortest 
one, and carry out the corresponding operation until we reach the 
end of the hydrodynamic time step t n+ . If in such a sub-cycle 
the shortest distance is d t , we translate the MCP by this distance 
Ti —> r; + dt n,, where ii; = Pi/pi is the propagation direction. 
We also attenuate its momentum according to Equation (38| l and ac¬ 
crue the momentum — Api, a and energy |Ap i>a |c transferred to the 
gas. If the shortest distance is d a , we do the same over this distance, 
but at the end of translation, we also randomize the MCP’s direc¬ 
tion n i —> n' and accrue the corresponding additional momentum 
—Api jS = pi (n' — n, ) and kinetic energy —v ■ Ap; jS transferred 
to gas. As a further variance-reduction tactic, given the statistical 
isotropy of n', we compute the momentum deposited in a scatter¬ 
ing event simply as —Api, s = — p,. If the shortest distance either 
is d a or db, we translate the MCP while attenuating its momentum 
and accruing the deposited energy and momentum. Then we either 
remove the MCP while instantaneously depositing the remaining 
momentum and energy to the gas (if the shortest distance is d a ), or 
transfer the MCP to its new host cell (or removed the MCP if it has 
reached a non-periodic boundary of the computational domain). 

As the MCPs are transported over a hydrodynamic time step 
At, the energy and momentum source terms for each cell are accu- 
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Figure 1. Spherically-averaged radiation energy density profiles in the 
three-dimensional radiative diffusion test (Section [3.1) . The analytical so¬ 
lutions (solid lines) and the numerical solutions (crosses) are shown at four 
different times, (0.2, 0.6, 1.2, 3.2) x 10 —10 s. The dashed line and the 
right axis show the refinement level of the AMR grid. 


mulated using 


cSo = AtAy ^(|A Pi , a |c v ■ Ap ilS ), 

(40) 

S - AtAU^ (APia + Api ’ s) ’ 

(41) 


where the sums are over all the absorption and scattering events that 
occurred in a specific computational cell during the hydrodynamic 
time step. The source terms are then substituted into Equations 0 
and ( |14[ > to compute the gas momentum and energy at the end of 
the hydrodynamic time step. 


3 ASSESSMENT OF NUMERICAL ALGORITHM 

To assess the validity of our radiation hydrodynamics implementa¬ 
tion, we performed a series of standard tests: a test of radiative dif¬ 
fusion in a scattering medium (Section [3.1| l, a test of gas-radiation 
thermal equilibration (Section[T2j, a test of thermal wave propaga¬ 
tion (Marshak wave; Section |3l3| l, and a radiative shock test (Sec¬ 
tion [T4]l. 

3.1 Radiative diffusion 

Here we test the spatial transport of MCPs across the AMR grid 
structure in the presence of scattering. In the optically thick limit, 
radiation transfer proceeds as a diffusion process. The setup is a cu¬ 
bical L — 1 cm 3 three-dimensional AMR grid with no absorption 
and a scattering coefficient of k s = 600 cm -1 . The scattering is as¬ 
sumed to be isotropic and elastic. We disable momentum exchange 
to preclude gas back-reaction and focus on testing the evolution of 
the radiation field on a non-uniform grid. 

At f = Os, we deposit an initial radiative energy (fjnit = 

3.2 x 10 6 erg) at the grid center in the form of 1,177,600 MCPs with 
isotropically sampled propagation directions. We lay an AMR grid 


Figure 2. Evolution of the gas (diamonds) and radiation (squares) energy 
density in the one-zone radiative equilibration test (Section [ v2| . The exact 
solution is shown with a solid (gas) and dashed (radiation) line. 


hierarchy such that the refinement level decreases with increasing 
radius as shown on the right axis of Figure [I] The grid spacing is 
Aa; = 2~ e ~ 2 L, where l is the local refinement level. A constant 
time step of At = 2x 10 -12 s is used and the simulation is run for 
4 x 10 _1 ° s. 

The diagnostic is the radiation energy density profile as a func¬ 
tion of distance from the grid center and time pe r ad(f, t). The exact 
solution in d spatial dimensions is given by 

where D = c/(dk s ) is the diffusion coefficient. 

Figure[T]shows the spherically-averaged radiation energy den¬ 
sity profile at four times. Excellent agreement of our MC results 
with the analytical expectation shows that our algorithm accurately 
captures radiation transport in a scattering medium. We have re¬ 
peated the test in one and two spatial dimensions and find the same 
excellent agreement. We have also checked that in multidimen¬ 
sional simulations, the radiation field as represented with MCPs 
preserves the initial rotational symmetry. 


3.2 Radiative equilibrium 


Here, in a one-zone setup, we test the radiation-gas thermal cou¬ 
pling in LTE. We enabled IMC with an im plic itness parameter of 
a = 1. Defining u r = aT 4 as in Section 


2.2 


where a is the ra¬ 
diation constant and T is the gas temperature, the stiff system of 
equations governing the gas and radiation internal energy density 
evolution is 


^ = fcaC (e ra d - , (43) 

fr-Kf - e, 4 1441 

where fc a is the absorption coefficient. We assume an ideal gas with 
adiabatic index 7 = 5/3. 
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Figure 3. Radiation energy density profiles in the Marshak wave test prob¬ 
lem at times 6 = (3, 10, 20) from left to right. Numerical integrations are 
shown by the data points. The solid lines are the reference solutions of |Su| 
|& Qlson|fT996) . 


We perform the one-zone test with parameters similar to those 
of |Tumer & Stone] ( |2001) and |Harries| ( (2011[ >, namely, the absorp¬ 
tion coefficient is fc a = 4.0 x 10 -s cm -i and the initial energy 
densities are pe = 10 8 erg cm” 8 and pe ra d = 0, respectively. The 
results and the corresponding exact solutions are shown in Figure 
[7] The MC solution agrees with the exact solutions within < 4% 
throughout the simulation. It shows that the physics of radiation- 
gas thermal exchange is captured well by our scheme, and that in 
static media, the scheme conserves energy exactly. 

As noted by |Cheatham| ( |2010) , the order of accuracy associ¬ 
ated with the IMC method depends on the choice of a and on the 
specifics of the model system. When the above test problem is re¬ 
peated with a = 0.5, the error is < 0.05% because the 0(At 2 )- 
residuals cancel out and the method is 0(At 3 ) -accurate. 


3.3 Marshak wave 

In this test, we simulate the propagation of a non-linear thermal 
wave, known as the Marshak wave, in a static medium in one spa¬ 
tial dimension jSu & Qlson|1996|[Gonzalez et al.|2007HKrumholz| 
|et al.|2 007; Zhang et ah 20 TT). The purpose of this standard test 
is to validate the code’s ability to treat nonlinear energy coupling 
between radiation and gas when the gas heat capacity is a function 
of gas temperature. We employ IMC with a = 1. 

Initially, a static, uniform slab with a temperature of 10 K oc¬ 
cupying the interval 0 ^ z ^ 15 cm is divided into 256 equal 
cells. An outflow boundary condition is used on the left and a re¬ 
flective one on the right. A constant incident flux Fine = ctsbTL 
of ksTinc = 1 keV thermal radiation, where asB and fee are the 
Stefan-Boltzmann and Boltzmann constants, is injected from the 
left (at z = 0). The gas is endowed with a constant, grey absorp¬ 
tion coefficient of fc a = 1 cm” 1 and a temperature-dependent vol¬ 
umetric heat capacity of c v = aT 3 . The constant a is related to the 
Su-Olson retardation parameter e via a = 4a/e and we set e = 1. 

The diagnostics for this test problem are the spatial radia¬ 
tion and gas energy density profiles at different times. [Su & Olson] 
( | 1996| > provided semi-analytical solutions in terms of the dimen¬ 
sionless position 

x = x/3 knZ (45) 


Figure 4. The same as Figure[3] but for the gas energy density. 


and time 

q _ 4ack a t 


(46) 


The radiation and gas internal energy density are expressed in terms 
of the dimensionless variables 


^ " inc 

c aT(x,8) 4 
[ ’ ] 4 Fine ’ 


(47) 

(48) 


where E r and T are the radiation energy density and gas temper¬ 
ature, respectively (note that the relation of v to the specific gas 
internal energy e is nonlinear). 

Figures [3] and [4] show the dimensionless numerical solution 
profiles at three different times, over-plotting the corresponding 
semi-analytical solutions of |Su & Qlson| fl996[ >. At early times, 
we observe relatively large deviations from the Su-Olson solutions 
near the thermal wavefront. This is in not surprising given that the 
solutions were obtained assuming pure radiative diffusion, yet at 
early times and near the wavefront, where the gas has not yet heated 
up, the optical depth is only about unity and the transport is not dif¬ 
fusive. |Gonzalez et al.| ( |2007] > observed the same early deviations 
in their computations based on the Ml closure. At later times when 
transport is diffusive, both the thermal wave propagation speed and 
the maximum energy density attained agree well with the Su-Olson 
solutions. 


3.4 Radiative shock 

To finally test the fully-coupled radiation hydrodynamics, we sim¬ 
ulate a radiative shock tube. As in the preceding tests, we use 
IMC with a — 1, but now, the gas is allowed to dynamically 
respond to the radiation. We adopt the setup and initial condi¬ 
tions of Ensman [1994| > and Commer gon et al.] < [201 1| > and simu¬ 
late both subcritical and supercritical shocks. The setup consists of 
a one-dimensional 7 x IO 10 cm-long domain containing an ideal 
gas with a mean molecular weight of p — 1 and adiabatic in¬ 
dex of 7 = 7/5. The domain is initialized with a uniform mass 
density po = 7.78 x 10” 1() gcm” 8 and a uniform temperature 
of To = 10 K. The gas has a constant absorption coefficient of 
k a = 3.1 x 10” 10 cm” 1 and a vanishing physical scattering coef¬ 
ficient. 

Initially, the gas is moving with a uniform velocity toward the 
left reflecting boundary. An outflow boundary condition is adopted 
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Table 1. Simulation parameters for the radiation-driven gaseous atmosphere. 


Run 

Initial Perturbation 

S 

(gem -2 ) 

9 

(10 —6 dyneg -1 ) 

/l* 

(10 -4 pc) 

u 

kyr 

T* 

/e,* 

^max /t* 

[Lx X Ly]/h* 

Ax/ht 

^max 

T10F0.02 

sin 

4.7 

37 

0.25 

0.045 

10 

0.02 

80 

512 x 256 

0.5 

7 

T03F0.50 

sin, x 

1.4 

1.5 

6.30 

1.1 

3 

0.5 

115 

512 x 2048 

1.0 

6 




Figure 5. Gas (solid curve) and radiation (dashed curve) temperature in the Figure 6. The same as Figure [5] but for the supercritical radiative shock 

subcritical radiative shock test at t = 3.8 X 10 4 s. The initial gas velocity test at t = 7.5 X 10 3 s and with initial velocity v = 20kms _1 . 

is v = 6 km s' and the profiles are plotted as a function of z = x — vt. 


on the right to allow inflow of gas at fixed density pa and fixed tem¬ 
perature To and also to allow the free escape of radiation MCPs. As 
gas collides with the reflecting boundary a shock wave starts propa¬ 
gating to the right. The thermal radiation in the compressed hot gas 
diffuses upstream and produces a warm radiative precursor. The 
shock becomes critical when the flux of thermal radiation is high 
enough to pre-heat the pre-shock gas to the post-shock temperature 
( jZerdovich & Raizer|1967| . We choose the incoming speed to be 
wo = 6 kms“ and 20 kms _1 in the subcritical and the supercriti¬ 
cal shock tests, respectively. 

|Mihalas & Mihala~s| l |1984} provide analytical estimates for the 
characteristic temperatures of the radiative shocks. For the subcrit¬ 
ical case, the post-shock temperature T 2 is estimated to be 


„ 2(7 - 1 )wq 

R (7 +l) 2 ' 


(49) 


Using the parameters for the subcritical setup, the analytical esti¬ 
mate gives T 2 ~ 812 K. In our simulation, the post-shock tem¬ 
perature at t = 3.8 x 10 4 s is T 2 — 800 K, which agrees with 
the analytical solution. The immediate pre-shock temperature T_ 
is estimated to be 


T- ~ 2( Z 1 } <t sb T 2 4 . (50) 

VoHpv 

Our simulation gives T- ~ 300 K while Tl is estimated to be 
T- = 270 K. Finally, the amplitude of the temperature spike can 


be estimated to be 

T+~T 2 + 5-^l-T_, (51) 

7+1 

which gives T+ ~ 990 K. It also close to the value we find, 
T + ~ 1000 K. In both cases, our simulations reproduce the ex¬ 
pected radiative precursors. Also, in the supercritical case, the pre¬ 
shock and the post-shock temperatures are identical, as expected. 


4 SETUP OF RADIATION-DRIVEN ATMOSPHERE 

We turn to the problem of how radiation drives an interstellar 
gaseous atmosphere in a vertical gravitational field. The problem 
was recently investigated by KT12 and KT13, by D14, and by 
RT15, using the FLD, VET, and Ml closure, respectively. Our aim 
is to attempt to reproduce these authors’ results, which are all based 
on low-order closures, using an independent method that does not 
rely on such a closure. Critical for the hydrodynamic impact of ra¬ 
diation pressure is the extent of the trapping of IR radiation by dusty 
gas. Therefore we specifically focus on the radiation transfer aspect 
of the problem and assume perfect thermal and dynamic coupling 
between gas and dust grains, T g = To = T and v g = = v. 

We follow the setup of KT12 and D14 as closely as possible. 
Taking that UV radiation from massive stars has been reprocessed 
into the IR at the source, we work in the grey approximation in 
which spectral averaging of the opacity is done only in the IR part 
of the spectrum. We set the Rosseland kr and Planck rep mean dust 
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opacities to 


« R ,p = (0.0316, 0.1) 


T 

10 A' 


2 

2 -1 
cm g 


(52) 


This model approximates a dusty gas in LTE at T ^ 150 K ( |Se-| 
|menov et al,|2003| l. Diverging slightly from KT12 and D14, who 
adopted the pure power-law scaling in Equation {52}, to approxi¬ 
mate the physical turnover in opacity, we cap both mean opacities 
to their values at 150 K above this threshold temperature. Overall, 
our opacity model is reasonable below the dust grain sublimation 
temperature ~ 1000 K. 

The simulation is set up on a two-dimensional Cartesian grid 
of size L x x L y . The grid is adaptively refined using the standard 
FLASH second derivative criterion in the gas density. The dusty 
gas is initialized as a stationary isothermal atmosphere. A time- 
independent, vertically incident radiation field is introduced at the 
base of the domain (y = 0) with a flux vector F* y. The gravita¬ 
tional acceleration is — gy. 

For notational convenience, we define a reference temperature 
T* = [Ft /(ca)] 1 / 4 , sound speed c* = scale 

height ht = c 2 /p, density p* = E//i* (where E is the ini¬ 
tial average gas surface density at the base of the domain), and 
sound crossing time t» = In the present setup F t = 

2.54 xlO 13 Lq kpc 2 and the mean molecular weight is p = 2.33 
as expected for molecular hydrogen with a 10 % helium molar frac¬ 
tion. The characteristic temperature is T* = 82 K and the corre¬ 
sponding Rosseland mean opacity is kr,* = 2.13 cm 2 g^ 1 . 

Following KT12 and KT13, we adopt two dimensionless pa¬ 
rameters to characterize the system: the Eddington ratio 


/e,» 


kr ,*F, 

gc 


(53) 


and the optical depth 


Tt = «R,*S. 


(54) 


The atmosphere is initialized at a uniform temperature Tt . The gas 
density is horizontally perturbed according to 


p(x,y) = 


x 


1 + 


1 + X 
4 


sin 



P te~ v/h % if e~ y/h * > W~ 10 , 
Pt 10 “ 10 , if e~ y ^ h * ^ 10 -1 °, 


(55) 


where \ x = 0.5 L x . D14 introduced x, a random variate uniformly 
distributed on [—0.25, 0.25], to provide an additional perturbation 
on top of the sinusoidal profile. If \ = 0, the initial density distri¬ 
bution reduces to that of KT12. 

KT12 found that at a given r*, the preceding setup has a hy¬ 
drostatic equilibrium solution when /e,* is below a certain critical 
value /E.crit- Note that KT12 defined /E,crit assuming the pure 
power-law opacity scaling in Equation {52}. With our capping of 
the opacities above 150 K, which breaks the dimensionless nature 
of the KT12 setup, the exact KT12 values for /E.crit cannot be 
directly transferred to our model. Nevertheless, we use their defini¬ 
tion of /E,crit simply to normalize our values of r* and /e,»- We 
attempt to reproduce the two runs performed by KT12. The first 
run T10F0.02 with t* = 10 and /e,* = 0.02 = 0.5 /E,crit lies in 
the regime in which such a hydrostatic equilibrium solution exists. 
The second run T03F0.5 with t« = 3 and /e,* = 0.5 = 3.8 /E,crit 
corresponds to the run performed by both KT12 and D14 that had 
the smallest ratio /e,* / /E.crit and was still unstable. The latter run 
probes the lower limit for the occurrence of a dynamically unstable 
coupling between radiation and gas. 



x/h. 


Figure 7. Gas density snapshots at four different times in the stable run 
T10F0.02. The fully simulation domain is larger than shown, 512x256 ti*; 
here, we only show the bottom quarter. The stable outcome of this run is 
consistent with the cited literature. 


The boundary conditions are periodic in the x direction, re¬ 
flecting at y = 0 (apart from the flux injection there), and out¬ 
flowing (vanishing perpendicular derivative) at y = L y . The re¬ 
flecting condition does not allow gas flow or escape of radiation. 
The outflow condition does allow free inflow or outflow of gas 
and escape of radiation. Unlike the cited treatments, we used non- 
uniform AMR. The AMR improves computational efficiency early 
in the simulation when dense gas occupies only a small portion of 
the simulated domain. In low-density cells, radiation streams al¬ 
most freely through gas; there, keeping mesh resolution low min¬ 
imizes the communication overhead associated with MCP han¬ 
dling while still preserving MCP kinematic accuracy. The appli¬ 
cation of AMR in conjunction with IMC is clearly not essential in 
two-dimensional, low-dynamic-range setups like the one presented 
here, but should become critical in three-dimensional simulations 
of massive star formation; thus, we are keen to begin validating it 
on simple test problems. 

To further economize computational resources, we require a 
density ^ 10 -6 p„ for thermal emission, absorption, and scattering 
calculations; below this density, the gas is assumed to be adiabatic 
and transparent. We also apply a temperature floor of 10 K. 

As the simulation proceeds, the total number of MCPs in¬ 
creases. To improve load balance, we limit the maximum number 
of MCPs allowed in a single computational block ( 8 x 8 cells) at 
the end of the time step to 64, or on average ~ 1 MCP per cell. (A 
much larger number of MCPs can traverse the block in the course 
of a time step.) If the number exceeds this specified maximum at 
the end of the radiation transport update, we merge some of the 
MCPs in a momentum- and energy-conserving fashion. We. how¬ 
ever, do not properly preserve spatial and higher-angular-moment 
statistical properties of the groups of MCPs subjected to merging. 
This deficiency is tolerable in the present simulation where merg¬ 
ing takes place only at the lowest level of refinement, where the 
radiation no longer affects the gas. In future applications, however, 
we will develop a manifestly more physical MCP merging strategy. 

The simulation parameters of the two runs are summarized 
in Table Q] The quoted cell width Ax is that at the highest level of 
mesh refinement (the cells are square). Gas with density > 10 -8 p« 
always resides at the maximum refinement level £ max throughout 
the simulation of duration f max . 
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Figure 9. Gas density in the unstable ran T03F0.50 at four different times as indicated. The panels display the full simulation domain. The flow morphology 
is consistent with that observed by D14 and RT15. Downward Rayleigh-Taylor plumes are visible in the second panel from the left. 


5 RESULTS 

5.1 Stable run T10F0.02 

Density snapshots at four different times in the simulation are 
shown in Figure]?] The simulation closely reproduces the quantita¬ 
tive results of both FLD (KT12) and VET (D14). Shortly after the 
beginning of the simulation, the trapping of radiation at the bottom 
of the domain by the dense dusty gas produces a rise in radiation 
energy density. As we assume LTE and perfect thermal coupling 
between gas and dust, the gas temperature increases accordingly. 
Specifically, after being heated up by the incoming radiation, the 
effective opacity at the midplane rises by a factor of ~ 10. The 
opacity rise enhances radiation trapping and the temperature rises 
still further to ~ (3 — 4) T, ~ 300 K. The heating drives the atmo¬ 
sphere to expand upward, but radiation pressure is not high enough 
to accelerate the slab against gravity. After the initial acceleration, 
the atmosphere deflates into an oscillatory, quasi-equilibrium state. 
This outcome is consistent with what has been found with FLD and 
VET. 


To better quantify how the dynamics and the degree of turbu¬ 
lence in the gas compare with the results of the preceding investi¬ 
gations, we compute the mass-weighted mean gas velocity 

1 r L y r L ->- 

(v) = — y y p{x, y)v(x, y)dxdy (56) 

and linear velocity dispersion 

1 fL y rLx 

ai=—J j p{x,y)(vi{x,y) - (vi)) dxdy, (57) 

where M is the total mass of the atmosphere and i indexes the 
coordinate direction. We also define the total velocity dispersion 

a = V a % +a l- 

The time evolution of (v y ) and the linear dispersions is shown 
in Figure[8] All the velocity moments are expressed as fractions of 
the initial isothermal sound speed c* = 0.54kms _1 . Both panels 
closely resemble those in Figure 2 of D14. Early on at ~10t*, ra¬ 
diation pressure accelerates the gas and drives growth in (v v ), a v , 
and a x . After this transient acceleration, ( v y ) executes dampled 
oscillations about zero velocity (the damping is likely of a numer¬ 
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Figure 8. Top panel: Time evolution of the mass-weighted mean velocity 
in the vertical direction in the stable run T10F0.02. Bottom panel: The cor¬ 
responding time evolution of the mass-weighted mean velocity dispersion. 
The linear dispersions a x and a y are shown with the dotted and dashed 
lines, respectively. 



t/t. 


Figure 10. Top panel: Time evolution of the mass-weighted mean velocity 
in the vertical direction in the unstable run T03F0.50 (black line). The col¬ 
ored lines are tracks from the cited references (see text and legend). Bottom 
panel: Mass-weighted mean velocity dispersions (see legend). In both pan¬ 
els, the late-time net acceleration and velocity dispersions are in agreement 
only with the results obtained by D14 with their short characteristics-based 
VET method. 

ical origin). The linear velocity dispersions also oscillate, but with 
smaller amplitudes < 0.4 c*. The oscillation period in cr^ is just 
slightly longer than that reported by D14. The agreement of our and 
VET results demonstrates the reliability of both radiation transfer 
methods. 

5.2 Unstable run T03F0.50 

To facilitate direct comparison with D14 and RT15, we introduce 
a random initial perturbation on top of the initial sinusoidal pertur¬ 
bation as in Equation ( |55[ >. The grid spacing Aat is twice of that 
adopted by KT12, D14, and RT15, but as we shall see, this coarser 
spacing is sufficient to reproduce the salient characteristics of the 
evolving system. MCP merging is activated at t = 36 f *. Figure 
[9] shows density snapshots at four different times. As in the stable 
case, the incoming radiation heats the gas at the bottom of the do¬ 



0 20 40 60 80 100 120 

i/t. 

Figure 11. Top panel: Time evolution of the volume-weighted Eddington 
ratio in the unstable run T03F0.50. The colored lines are tracks from the 
cited references (see text and legend). Middle panel: The volume-weighted 
mean total vertical optical depth. Bottom panel: The flux-weighted mean 
optical depth. 

main and opacity jumps. The flux soon becomes super-Eddington 
and a slab of gas is lifted upward. At 39t*, fragmentation of the 
slab by the RTI is apparent; most the gas mass becomes concen¬ 
trated in dense clumps. 

We note that the slab lifting and the subsequent fragmenta¬ 
tion are consistently observed in all radiative transfer approaches; 
differences become apparent only in the long-term evolution. As 
in the VET simulation (D14), coherent gaseous structures in our 
simulation continue to be disrupted and accelerated. Qualitatively, 
radiation drives gas into dense, low-filling-factor filaments embed¬ 
ded in low density (10 -3 — 10 -4 ) p* gas. At 115f*, the bulk of 
the gas has a net upward velocity and has been raised to altitudes 
y ~ 1500/1*. 

Figure [TO] compares the time evolution of the bulk velocity 
(v y ) and velocity dispersions a x ^ y with the corresponding tracks 
from the published FLD/VET and Ml simulations (respectively, 
D14 and RT15). Initially, (v y ) rises steeply as the gas slab heats up 
and the incoming flux becomes super-Eddington. All simulations 
except for the one performed with the Ml closure without radia¬ 
tion trapping (RTI5) exhibit a similar initial rise. At ~ 25 i*, the 
RTI sets in and the resulting filamentation reduces the degree of 
radiation trapping. This in turn leads to a drop in radiation pres¬ 
sure and (v y ) damps down under gravity. The transient rise and 
drop in (v y ) is observed with all the radiative transfer methods, 
although the specific times of the acceleration-to-deceleration tran¬ 
sition differ slightly. The bulk velocity peaks at ( v y } ~ 12 c* in 
IMC and at ~ 9 c* in VET. The subsequent kinematics differs sig¬ 
nificantly between the methods. In IMC and VET, the gas filaments 
rearrange in a way that enables resumption of upward acceleration 
after ~ (50 — 60) t*. At late times, the secondary rise in (v y ) does 
not seem to saturate in IMC as it does in VET. Otherwise, the IMC 
and VET tracks are very similar to each other. In FLD and Ml, 
however, gas is not re-accelerated after the initial transient accel¬ 
eration. Instead, it reaches a turbulent quasi-steady state in which 
gas is gravitationally confined at the bottom of the domain and (v y ) 
fluctuates around zero. 

The evolution of velocity dispersions in IMC is also in close 
agreement with VET. Before the RTI onset, the dispersions rise 
slightly to a y > 1 c*. Once the RTI develops and the slab frag¬ 
ments, a y increases rapidly and u x somewhat more gradually. A 
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drop in a y is observed at ~ 75 t», but after that time, the vertical 
dispersion rises without hints of saturation. Velocity dispersions at 
the end of our simulations are consistent with those in VET. In FLD 
and Ml, on the other hand, the asymptotic turbulent quasi-steady 
states have smaller velocity dispersions ~ 5 c*. 

To further investigate the coupling of gas and radiation, we 
follow KT12 and KT13 to define three volume-weighted quantities: 
the Eddington ratio 


/e,v = 


{K R pFy)y 

cgp 


the mean total vertical optical depth 

Tv = L v (krp)v, 

and the flux-weighted mean optical depth 

(tvR pFy) V 


TP = Ly 


{Fyh 


(58) 


(59) 


(60) 


•>v = 


where F y is the flux in the y direction and 
L” 1 Ly 1 j 0 v f x • dxdy denotes volume avearge. 

Figure ED compares the time evolution of the volume- 
weighted quantities in our IMC run with those in FLD, Ml, and 
VET. The evolution of /e,v in IMC matches both qualitatively and 
quantitatively that in VET over the entire course of the run. Com¬ 
mon to all the simulations except the one carried out with the Ml 
closure without radiative trapping, the mean Eddington ratio in¬ 
creases from its initial value of /e,v = 0.5 to super-Eddington 
values soon after the simulation beginning. Then it immediately 
declines toward /e,v 1-5- After ~ 201*, all methods become 
sub-Eddington, with the Ml with radiative trapping and the FLD 
exhibiting the most significant decline. Then beyond ~ 60 f*, all 
simulations attain near-unity Eddington ratios. D14 pointed out that 
the time evolution of (v y ) is in general sensitive to the value of 
/e, v, namely, increases when J' e, v 1 and decreases other 
wise. It is observed that IMC stays slightly super-Eddington at late 
times, similar to VET. The observed continuous acceleration of the 
gas with IMC suggests that gas dynamics can be very different be¬ 
tween simulations with similar volume-average Eddington rations 
as long as the simulations are performed with different radiative 
transfer methods. 

The middle panel of Figure im shows the evolution of the 
volume-weighted mean total vertical optical depth tv- Since this 
quantity depends only on the gas state but not on the noisier radi¬ 
ation state, the IMC track is smooth. It is a global estimate of the 
optical thickness of the gas layer and we expect its behavior to be 
related to that of /e,v- The IMC track seems a flattened and down- 
scaled version of the others. This should be an artifact of the precise 
choice of the opacity law. We cap the opacities kr,p at their val- 


6 CONCLUSIONS 

We applied the Implicit Monte Carlo radiative transfer method to a 
standard two-dimensional test problem modeling the radiation hy¬ 
drodynamics of a dusty atmosphere that is accelerated against grav¬ 
ity by an IR radiation field. The atmosphere is marginally capable 
of trapping the transiting radiation. We consider this idealized sim¬ 
ulation a necessary stepping stone toward characterizing the dy¬ 
namical impact of the radiation emitted by massive stars and active 
galactic nuclei. We compare our IMC-derived results with those us¬ 
ing low-order closures of the radiative transfer hierarchy that have 
been published by other groups. Our particle-based approach en¬ 
ables independent validation of the hitherto tested methods. 

Sufficiently strong radiation fluxes universally render the at¬ 
mosphere turbulent, but its bulk kinematics differs between the 
VET and IMC methods on the one hand and the FLD and Ml meth¬ 
ods on the other. We find that the former continue to accelerate the 
atmosphere against gravity in the same setup in which the latter reg¬ 
ulate the atmosphere into a gravitationally-confined, quasi-steady 
state. This exposes shortcomings of the local closures. Namely, in 
complex geometries, the FLD seems to allow the radiation to more 
easily escape through optically thin channels. This can be under¬ 
stood in terms of a de facto artificial re-collimation of the radiation 
field diffusing into narrow, optically-thin channels from their more 
optically thick channel walls. In the limit in which the radiation 
freely streams in the channels, the flux in the channels becomes 
equal to what it would be for a radiation field in which the photon 
momenta are aligned with the channel direction. Indeed, D14 argue 
that in the optically thin regime, the FLD's construction of the radi¬ 
ation flux is inaccurate in both its magnitude and direction, and has 
the tendency to reinforce the formation of such radiation-leaking 
channels. 

Whether outflowing or gravitationally-confined, the turbulent 
atmosphere seems to reach a state approximately saturating the Ed¬ 
dington limit. The nonlinearity arising from the increase of dust 
opacity with temperature introduces the potential for bi-stability in 
the global configuration. Subtle differences between numerical clo¬ 
sures can be sufficient to force the solution into degenerate, qual¬ 
itatively different configurations. Robust radiation-hydrodynamic 
modeling seems to demand redundant treatment with distinct nu¬ 
merical methods including the IMC. 

Future work will of course turn to more realistic astrophysical 
systems. For example, the role of radiation trapping and pressure in 
massive star forming regions remains a key open problem, both in 
the context of the nearby jKrumholz & Matzner|2009 Krumholz 


|et al.||2012| [2014] |Coker etah 


2013} jLopez et al.|[2014| and the 


distant i Riechers & et al.||2013 > universe. Radiative reprocessing 


by photoionization and dust requires a frequency-resolved treat¬ 
ment of the radiation field as well as a generalization the IMC 
method to nonthermal processes. The assumption of perfect gas- 
dust thermal coupling can be invalid and the respective temper¬ 
atures must be hacked separately. Numerical treatments may be 


ues at T = 150 K, whereas the other authors allow the kckT 

required to resolve dust sublimation fronts ([Kuiper et al. 

2010 

scaling to extend at T > 150 K. Therefore, our choice of opac- 

and radiation pressure on metal lines {Tanaka & Nakamoto 

2011 


ity underestimates the strength of radiation pressure compared to 
the cited studies, but this discepancy does not appear to affect the 
hydrodynamic response of the gas. 


massive- 


The bottom panel of Figure 11 shows the ratio tf/tv- Note 


that tf is the true effective optical depth felt by the radiation. There¬ 
fore, a small tf/tv implies a higher degree of flux-density anti¬ 
correlation. The evolution of this ratio is similar in all radiation 
transfer methods. 


star-forming cores, multifrequency radiative transfer may be of 
essence for robust estimation of the final characteristic stellar mass 
scale and the astronomically measurable accretion rate ([Yorke &| 
|Sonnhalter||2002[ |Tan et al.|2014| >. Photoionization can set the fi¬ 
nal stellar masses through fragmentation-induced starvation E 
t ers et al.|2010| l. The star formation phenomenon spans a huge dy¬ 
namic range that can be effectively treated with telescopic AMR 
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grids constructed to ensure that the local Jeans length is always ad¬ 
equately resolved. It will likely be necessary to invent new acceler¬ 
ation schemes for improving the IMC method’s efficiency in such 
heterogeneous environments. One promising direction is the intro¬ 
duction of MCP splitting (see, e.g., |Harries|2015[ where MCP split¬ 
ting is applied in methods developed to simulate radiation transfer 
in massive star forming systems). 
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